%--------------------------------------------------------------------------
% function computes the stabilization times numerically and using the
% multiple scales approach
%--------------------------------------------------------------------------


function [t_num, t_asy_max] = comp_times(k, p)

if nargin == 1
    p = params;
end

% the numeric stabilization time
s = stab(k, p);
A = @(t) interp1(s.x, max(abs(s.y)), t);
t_num = fminsearch(@(t) -A(t), 2);
fprintf('t_num = %.4e\n', t_num);


% time predicted by neutral stab curves crossing
t_ns = fzero(@(t) p.Ma - comp_Ma(t, k), t_num);
fprintf('t_ns = %.4e\n', t_ns);

% time predicted by eigenvalue
h = -(1 + lambertw(-p.beta / (p.beta - 1) * exp(((-p.beta + p.delta * s.x) / (p.beta - 1))))) * (p.beta - 1);
for i = 1:length(s.x)
    ev(i) = comp_eigs(k, p, h(i));
end
f = @(t) interp1(s.x, ev, t);
t_zero_eig = fzero(@(t) f(t), t_num);
fprintf('t_zero_eig = %.4e\n', t_zero_eig);

% time predicted from asymptotic amplitude
sigma = @(t) interp1(s.x, ev, t);
A_asy = @(T) exp(quad(@(t) sigma(t), 0, T));
t_asy_max = fminsearch(@(t) -A_asy(t), t_num);
fprintf('t_asy_max = %.4e\n', t_asy_max);

function Ma = comp_Ma(t, k)

[k, Ma] = exact_neutral_stab(k, t, 0);